###############################################################################################
# 
# Sensitivity Analysis: Unobserved Confounding
# 
###############################################################################################

#########################################################################################
### Democratic Republic of Congo

est <- m.out.congo$par.outcome[14]      # Estimate of Treatment Effect
se <- m.out.congo$se.outcome[14]           # SE of Treatment Effect
df <- 919 - 43 - 1     # Degrees of Freedom

bias <- function(R2.yz, R2.dz){  
  bias = se * sqrt( (R2.yz * R2.dz)/(1-R2.dz) * df )
  return(est-bias)
}

x <- seq(0, 1, by=.01)
y <- seq(0, 1, by=.01)
sens <- outer(x, y, FUN=bias)

lm.drc <- ictreg(Y ~ female + age.z + edu.z + income.z  + hh_size.z +  murder_yes 
               + terr_2 + terr_3 + terr_4 + terr_5 + terr_6 + activity_prev, treat="treatment", J=3, data=congo, method="lm")

t <- lm.drc$par.treat/lm.drc$se.treat
p.R2.dz <- t^2/(t^2 + lm.drc$resid.df)

t <- m.out.congo$par.outcome/m.out.congo$se.outcome
p.R2.yz <- t^2/(t^2 + df)
p.R2.yz <- p.R2.yz[-c(length(p.R2.yz))]

par(mfrow=c(3,1), mar=c(4, 3, 2, 2), cex.lab=.8, cex.axis=.8, las=1, mgp=c(1.7, .5, 0))
contour(sens, nlevels=50, xlab=expression(paste("Hypothetical Partial ", R^2, " of Unobserved Confounder(s) with Sexual Violence")), 
        ylab=expression(paste("Hypothetical Partial ", R^2, "of Unobserved Confounder(s) with Participation")),
        labcex=1, method="edge", axes=F, col="red", lwd=1, xlim=c(0, .2), ylim=c(0, .25), main="DRC")
#contour(sens, levels=0, labcex=1, col="red", lwd=2, add=T)
axis(1, col="white")
axis(2, col="white")
#grid(col="grey", lwd=.5, lty=1)
box(col="darkgrey")

points(p.R2.dz[-1], p.R2.yz[-1], pch=19, cex=.5)

points(p.R2.dz[13], p.R2.yz[13], pch=19, cex=.5, col="orange")
text(p.R2.dz[13], p.R2.yz[13], "Pre-exposure Participation", cex=.6, col="orange", pos=4)

points(p.R2.dz[13]*2, p.R2.yz[13]*2, pch=19, cex=.5, col="orange")
text(p.R2.dz[13]*2, p.R2.yz[13]*2, "Pre-exposure Participation x 2", cex=.6, col="orange", pos=4)

points(p.R2.dz[13]*4, p.R2.yz[13]*4, pch=19, cex=.5, col="orange")
text(p.R2.dz[13]*4, p.R2.yz[13]*4, "Pre-exposure Participation x 4", cex=.6, col="orange", pos=4)

points(p.R2.dz[4], p.R2.yz[4], pch=19, cex=.5, col="purple")
points(p.R2.dz[4]*5, p.R2.yz[4]*5, pch=19, cex=.5, col="purple")
points(p.R2.dz[4]*10, p.R2.yz[4]*10, pch=19, cex=.5, col="purple")
text(p.R2.dz[4], p.R2.yz[4], "Education", cex=.6, col="purple", pos=4)
text(p.R2.dz[4]*5, p.R2.yz[4]*5, "Education x 5", cex=.6, col="purple", pos=4)
text(p.R2.dz[4]*10, p.R2.yz[4]*10, "Education x 10",  cex=.6, col="purple", pos=4)

text(.12, .07, "Effect would be zero.",  cex=.6, col="red", pos=4)
arrows(.12, .07, .1, .05, length=.03, col="red")

#########################################################################################
### Liberia

est <- m.out.liberia$par.outcome[10]      # Estimate of Treatment Effect
se <- m.out.liberia$se.outcome[10]           # SE of Treatment Effect
df <- 6361 - 31 - 1     # Degrees of Freedom

bias <- function(R2.yz, R2.dz){  
  bias = se * sqrt( (R2.yz * R2.dz)/(1-R2.dz) * df )
  return(est-bias)
}

x <- seq(0, 1, by=.01)
y <- seq(0, 1, by=.01)
sens <- outer(x, y, FUN=bias)

lm.lib <- ictreg(Y ~ female + age.z + edu.z + income.z + hh_size.z + cw_kill +  as.factor(county), treat="treatment", J=3, data=liberia, method="lm")

t <- lm.lib$par.treat/lm.lib$se.treat
p.R2.dz <- t^2/(t^2 + lm.lib$resid.df)

t <- m.out.liberia$par.outcome/m.out.liberia$se.outcome
p.R2.yz <- t^2/(t^2 + df)
p.R2.yz <- p.R2.yz[-c(length(p.R2.yz))]

par(mar=c(4, 3, 2, 2), cex.lab=.8, cex.axis=.8, las=1, mgp=c(1.7, .5, 0))
contour(sens, nlevels=50, xlab=expression(paste("Hypothetical Partial ", R^2, " of Unobserved Confounder(s) with Sexual Violence")), 
        ylab=expression(paste("Hypothetical Partial ", R^2, "of Unobserved Confounder(s) with Participation")),
        labcex=1, method="edge", axes=F, col="red", lwd=1, xlim=c(0, .2), ylim=c(0, .25), main="Liberia")
#contour(sens, levels=0, labcex=1, col="red", lwd=2, add=T)
axis(1, col="white")
axis(2, col="white")
#grid(col="grey", lwd=.5, lty=1)
box(col="darkgrey")

points(p.R2.dz[-1], p.R2.yz[-1], pch=19, cex=.5)

points(p.R2.dz[4], p.R2.yz[4], pch=19, cex=.5, col="purple")
points(p.R2.dz[4]*15, p.R2.yz[4]*15, pch=19, cex=.5, col="purple")
points(p.R2.dz[4]*30, p.R2.yz[4]*30, pch=19, cex=.5, col="purple")

text(p.R2.dz[4], p.R2.yz[4], "Education", cex=.6, col="purple", pos=4)
text(p.R2.dz[4]*15, p.R2.yz[4]*15, "Education x 15", cex=.6, col="purple", pos=4)
text(p.R2.dz[4]*30, p.R2.yz[4]*30, "Education x 30",  cex=.6, col="purple", pos=4)

text(.12, .03, "Effect would be zero.",  cex=.6, col="red", pos=4)
arrows(.12, .03, .1, .01, length=.03, col="red")

#########################################################################################
### Sri Lanka

est <- m.out.sri$par.outcome[18]      # Estimate of Treatment Effect
se <- m.out.sri$se.outcome[18]           # SE of Treatment Effect
df <- 920 - 52 - 1     # Degrees of Freedom

bias <- function(R2.yz, R2.dz){  
  bias = se * sqrt( (R2.yz * R2.dz)/(1-R2.dz) * df )
  return(est-bias)
}

x <- seq(0, 1, by=.01)
y <- seq(0, 1, by=.01)
sens <- outer(x, y, FUN=bias)

lm.sri <- ictreg(Y ~ female + age.z + edu.z + income.z + hh_size.z +  killed + trauma + prov.2 + prov.3 + prov.4 + prov.5 + prov.6 + prov.7 + prov.8 + prov.9 + prior, treat="treatment", J=3, data=sri)

t <- lm.sri$par.treat/lm.sri$se.treat
p.R2.dz <- t^2/(t^2 + 884)

t <- m.out.sri$par.outcome/m.out.sri$se.outcome
p.R2.yz <- t^2/(t^2 + df)
p.R2.yz <- p.R2.yz[-c(length(p.R2.yz))]

par(mar=c(4, 3, 2, 2), cex.lab=.8, cex.axis=.8, las=1, mgp=c(1.7, .5, 0))
contour(sens, nlevels=50, xlab=expression(paste("Hypothetical Partial ", R^2, " of Unobserved Confounder(s) with Sexual Violence")), 
        ylab=expression(paste("Hypothetical Partial ", R^2, "of Unobserved Confounder(s) with Participation")),
        labcex=1, method="edge", axes=F, col="red", lwd=1, xlim=c(0, .2), ylim=c(0, .25), main="Sri Lanka")
#contour(sens, levels=0, labcex=1, col="red", lwd=2, add=T)
axis(1, col="white")
axis(2, col="white")
#grid(col="grey", lwd=.5, lty=1)
box(col="darkgrey")

points(p.R2.dz[-1], p.R2.yz[-1], pch=19, cex=.5)

points(p.R2.dz[4], p.R2.yz[4], pch=19, cex=.5, col="purple")
points(p.R2.dz[4]*5, p.R2.yz[4]*5, pch=19, cex=.5, col="purple")
points(p.R2.dz[4]*10, p.R2.yz[4]*10, pch=19, cex=.5, col="purple")

text(p.R2.dz[4], p.R2.yz[4], "Education", cex=.6, col="purple", pos=4)
text(p.R2.dz[4]*5, p.R2.yz[4]*5, "Education x 5", cex=.6, col="purple", pos=4)
text(p.R2.dz[4]*10, p.R2.yz[4]*10, "Education x 10",  cex=.6, col="purple", pos=4)

points(p.R2.dz[17], p.R2.yz[17], pch=19, cex=.5, col="orange")
points(p.R2.dz[17]*15, p.R2.yz[17]*15, pch=19, cex=.5, col="orange")
points(p.R2.dz[17]*30, p.R2.yz[17]*30, pch=19, cex=.5, col="orange")
points(p.R2.dz[17]*45, p.R2.yz[17]*45, pch=19, cex=.5, col="orange")

text(p.R2.dz[17], p.R2.yz[17], "Pre-exposure  Participation",  cex=.6, col="orange", pos=4)
text(p.R2.dz[17]*15, p.R2.yz[17]*15, "Pre-exposure  Participation x 15",  cex=.6, col="orange", pos=4)
text(p.R2.dz[17]*30, p.R2.yz[17]*30, "Pre-exposure  Participation x 30",  cex=.6, col="orange", pos=4)
text(p.R2.dz[17]*45, p.R2.yz[17]*45, "Pre-exposure  Participation x 45",  cex=.6, col="orange", pos=4)

text(.1, .13, "Effect would be zero.",  cex=.6, col="red", pos=4)
arrows(.1, .13, .08, .11, length=.03, col="red")

#########################################################################################

